library(haven)
library(ggplot2)
library(tidyverse)
library(patchwork)
library(kableExtra)
library(stargazer)
library(lmtest)
library(plm)
library(caschrono)
library(astsa)
library(forecast)
library(DT)
library(knitr)
library(plotly)
library(dplyr)
data <- read_dta("data_V3.dta")
data <- data[data$age != 996 & data$educyr <= 90,]
Le but de cet exercice est d’étudier les déterminants des dépenses en soins dentaires. Pour ce faire, nous disposons d’une base de données de panel. Elle a été construite à partir des réponses de 73 017 individus vivant aux Etats-Unis et qui ont été interrogés en 2015 et 2019.
Les variables sont les suivantes :
Avant de commencer notre travail, nous avons étudié longuement la base de données afin de supprimer les variables aberrantes telles que des individus qui avaient plusieurs centaines d’années (ou bien certains qui avaient changé de sexe). Ce sont, a priori, des erreurs de saisie. Nous avons également modifié certaines variables afin de faciliter notre étude.
Il est possible d’analyser les déterminants des dépenses en soins dentaires pour déterminer s’il existe des relations avec divers facteurs tels que le revenu, le sexe, le niveau d’éducation, l’âge, la possession d’une assurance maladie et la présence d’une assurance Medicare. Cette analyse peut fournir des informations sur les inégalités dans l’accès aux soins dentaires. Elle peut également aider à comprendre comment les différences socio-économiques influencent les dépenses en soins dentaires.
Les individus dans notre base de données ont en moyenne 48 ans et sont tous des adultes. Ils ont généralement passé entre 12 et 16 années à étudier. Leurs revenus sont très variés, ce qui nous indique que les individus proviennent de plusieurs classes sociales différentes. En moyenne, les dépenses en soins dentaires sont de 307,5 dollars par an et par personne, mais cette moyenne est influencée par des valeurs extrêmement élevées. La médiane de ces dépenses est de 0 dollars, ce qui semble cohérent car peu de personne ont des dépenses dentaires, mais lorsqu’une personne doit payer pour des soins dentaires, cela coûte cher.
summary(data[,c("age","educyr","inctot","dvexptot")]) %>%
kbl() %>%
kable_minimal()
| age | educyr | inctot | dvexptot | |
|---|---|---|---|---|
| Min. :19.00 | Min. : 0.00 | Min. :-309948 | Min. : 0.0 | |
| 1st Qu.:33.00 | 1st Qu.:12.00 | 1st Qu.: 10800 | 1st Qu.: 0.0 | |
| Median :48.00 | Median :13.00 | Median : 25021 | Median : 0.0 | |
| Mean :48.38 | Mean :13.02 | Mean : 36915 | Mean : 307.5 | |
| 3rd Qu.:62.00 | 3rd Qu.:16.00 | 3rd Qu.: 50000 | 3rd Qu.: 200.0 | |
| Max. :85.00 | Max. :17.00 | Max. : 731653 | Max. :81000.0 |
Il y a 54% de femmes et 46% d’hommes dans la base de données. Il y a donc un peu plus de femmes que d’hommes, mais le déséquilibre n’est pas très important. Selon les données de la Banque Mondiale, la part de femme aux USA varie entre 50 et 51% entre 1960 et 2021. (source).
Un quart des individus présent dans la base de données ont l’assurance medicare, c’est un programme de santé publique visant les personnes handicapées et les personnes âgées afin de leur donner une couverture medicale.
Nous avons 11% des individus qui n’ont pas d’assurance santé.
(ggplot(data, aes(x = sex)) +
geom_bar(color = "black", fill = "#267334", alpha = 0.8, width=0.6) +
ggtitle("Sexe") +
labs(subtitle = "1 : homme et 2 : femme")+
ylab("Effectifs")+
theme_bw()+
ggplot(data, aes(x = year)) +
geom_bar(color = "black", fill = "#267334", alpha = 0.8) +
ggtitle("Année") +
ylab("Effectifs")+
theme_bw())/
(ggplot(data, aes(x = himcare)) +
geom_bar(color = "black", fill = "#267334", alpha = 0.8, width=0.6) +
ggtitle("Assurance medicare") +
labs(subtitle = "1 non et 2 oui")+
ylab("Effectifs")+
theme_bw()+
ggplot(data, aes(x = hinotcov)) +
geom_bar(color = "black", fill = "#267334", alpha = 0.8, width=0.6) +
ggtitle("N’a pas d’assurance santé") +
labs(subtitle = "1 non et 2 oui")+
ylab("Effectifs")+
theme_bw())
Nous pouvons constater que 89% des individus présents dans la base de données ont déjà travaillé. On a aussi la surprise de voir qu’une petite poignée d’individus ne savent pas s’ils ont déjà travailler ou non.
tab <- as.matrix(table(data[,c("workev")]))
rownames(tab) <- c("non","oui","a refusé de répondre","non vérifié","ne sait pas")
tab <- cbind(tab,tab/sum(tab)*100)
colnames(tab) <- c("Effectifs","Proportions (%)")
tab %>% round(2)%>%
kbl(caption = "A déjà travaillé ?", booktabs = T, linesep = "") %>%
kable_minimal(full_width = F)
| Effectifs | Proportions (%) | |
|---|---|---|
| non | 11018 | 10.60 |
| oui | 92102 | 88.58 |
| a refusé de répondre | 568 | 0.55 |
| non vérifié | 25 | 0.02 |
| ne sait pas | 266 | 0.26 |
Par la suite, nous allons modifier la varible workev,
nous garderons seulement les réponses oui et
non, elle fonctionnera donc comme les autres
indicatrices.
\(dvexptot_{it}= \beta_0 + \beta_1age_{it} + \beta_2sex_{it} + \beta_3educyr_{i} + \beta_4inctot_{it} + \beta_5hinotcov_{it} + \beta_6workev_i + \beta_7himcare_{it} + \upsilon_{it}\)
avec \(\upsilon_{it} = \epsilon_{it} + c_i\)
data$sex <- as.factor(data$sex)
data$himcare <- as.factor(data$himcare)
data$hinotcov <- as.factor(data$hinotcov)
data$workev <- as.factor(data$workev)
library(stats)
data <- data[data$age != 996 & data$educyr <= 90 & (data$workev == 1 | data$workev == 2),]
pool <- plm(dvexptot ~ age + sex + educyr + inctot + hinotcov + workev + himcare, data = data,
model="pooling", index=c("id","year"))
pool_cor <- coeftest(pool, vcov=vcovHC(pool, method="arellano"))
En corrigeant les écarts-types, nous avons amélioré la
précision de nos estimations, à l’exception de la variable
himcare.
De plus, nous constatons que presque tous les coefficients du modèle
sont significatifs au seuil de 1%, ce qui indique que
chaque variable explicative a un impact statistiquement significatif sur
la variable dépendante. workev est la seule variable non
significative.
| Variable dépendante : dvexptot | ||
| Modèle poolé | Modèle poolé corrigé | |
| age | 4.159*** | 4.159*** |
| (0.283) | (0.283) | |
| sex2 | 63.037*** | 63.037*** |
| (7.089) | (7.492) | |
| educyr | 26.216*** | 26.216*** |
| (1.280) | (1.264) | |
| inctot | 0.002*** | 0.002*** |
| (0.0001) | (0.0001) | |
| hinotcov2 | -97.688*** | -97.688*** |
| (11.664) | (7.426) | |
| workev2 | 6.474 | 6.474 |
| (11.752) | (11.990) | |
| himcare2 | 72.639*** | 72.639*** |
| (11.882) | (14.296) | |
| Constant | -349.614*** | -349.614*** |
| (22.216) | (22.070) | |
| Observations | 103,120 | |
| R2 | 0.024 | |
| R2 ajusté | 0.024 | |
| F Statistic | 363.255*** (df = 7; 103112) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
Notre individu de référence est un homme qui n’a pas
d’assurance Medicare, qui a une assurance santé privée et qui n’a jamais
travaillé. Le coefficient associé à himcare est positif,
montrant que l’assurance Medicare peut augmenter les dépenses de
santé dentaire. Le coefficient de hinotcov est
négatif, reflétant l’impact négatif de ne pas avoir d’assurance
santé sur les dépenses de santé. Les années d’études, l’âge, un
revenu élevé et l’expérience professionnelle sont des facteurs qui
augmentent les dépenses de santé. Etre une femme a également un impact
positif sur les dépenses de santé.
\(dvexptot_{it}= \beta_0 + \beta_1age_{it} + \beta_2sex_{it} + \beta_3educyr_{i} + \beta_4inctot_{it} +\beta_5hinotcov_{it} + \beta_6workev_i + \beta_7himcare_{it} + c_i + \epsilon_{it}\)
pool_ea <- plm(dvexptot ~ age + sex + educyr + inctot + hinotcov + workev + himcare, data = data,
model="random", effect="individual",index=c("id","year"))
pool_ea_corr <- coeftest(pool_ea, vcov=vcovHC(pool_ea, method="arellano"))
| Variable dépendante : dvexptot | ||
| Modèle EA | Modèle EA corrigé | |
| age | 4.104*** | 4.104*** |
| (0.297) | (0.284) | |
| sex2 | 62.651*** | 62.651*** |
| (7.462) | (7.476) | |
| educyr | 26.313*** | 26.313*** |
| (1.345) | (1.262) | |
| inctot | 0.002*** | 0.002*** |
| (0.0001) | (0.0001) | |
| hinotcov2 | -97.005*** | -97.005*** |
| (12.096) | (7.553) | |
| workev2 | 6.564 | 6.564 |
| (12.315) | (12.008) | |
| himcare2 | 74.702*** | 74.702*** |
| (12.452) | (14.360) | |
| Constant | -347.061*** | -347.061*** |
| (23.305) | (22.091) | |
| Observations | 103,120 | |
| R2 | 0.022 | |
| R2 ajusté | 0.022 | |
| F Statistic | 2,284.897*** | |
| Note: | p<0.1; p<0.05; p<0.01 | |
| Modèle poolé | Modèle EA | |
| age | 4.159*** | 4.104*** |
| (0.283) | (0.284) | |
| sex2 | 63.037*** | 62.651*** |
| (7.492) | (7.476) | |
| educyr | 26.216*** | 26.313*** |
| (1.264) | (1.262) | |
| inctot | 0.002*** | 0.002*** |
| (0.0001) | (0.0001) | |
| hinotcov2 | -97.688*** | -97.005*** |
| (7.426) | (7.553) | |
| workev2 | 6.474 | 6.564 |
| (11.990) | (12.008) | |
| himcare2 | 72.639*** | 74.702*** |
| (14.296) | (14.360) | |
| Constant | -349.614*** | -347.061*** |
| (22.070) | (22.091) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
Il y a peu de différence entre le modèle poolé et le modèle à effets aléatoires. Les signes, les écarts-types et les coefficients sont similaires dans les deux modèles. Cependant, le choix peut se porter sur un modèle ou sur l’autre en fonction des hypothèses de base. En effet, un modèle poolé ou un modèle à effets aléatoires doit respecter la condition \(Cov(c_i,x_it) = 0\). Mais on sait que l’estimateur des MCO sur le modèle poolé est convergent (si la condition précédente est vérifiée) mais non efficace, c’est pourquoi un modèle à effets aléatoires est préférable.
Notons que l’on peut considérer que le revenu total d’un individu peut être corrélé avec le nombre d’années d’éducation qu’il a reçues. Dans ce cas, on est un peu embêté pour choisir entre ces deux modèles car ils sont probablement non convergents. Cependant dans le cas d’un modèle à effets fixes nous pouvons avoir un modèle convergent dans le cas où \(Cov(c_i,x_it) ≠ 0\)
Nous allons effectuer un modèle within. Comme nous sommes dans un modèle à effets fixes on s’attend à ce que nous n’ayons pas de coefficients pour les variables qui ne varient pas dans le temps.
\(dvexptot_{it} - \overline{dvexptot_{it}} = \beta_1(age_{it} - \overline{age_{it}}) + \beta_2(sex_{it} - \overline{sex_{it}}) + \beta_3(educyr_{i} - \overline{educyr_{i}})+ \beta_4(inctot_{it} - \overline{inctot_{it}})\) \(+ \beta_5(hinotcov_{it}- \overline{hinotcov_{it}}) + \beta_6(workev_i - \overline{workev_it}) + \beta_7(himcare_{it} - \overline{workev_i}) + (\epsilon_{it} - \overline{\epsilon_{it}})\)
fixed <- plm(dvexptot ~ age + sex + educyr + inctot + hinotcov + workev + himcare, data = data,
model="within", index = c("id", "year"))
fixed_cor <- coeftest(fixed, vcov=vcovHC(fixed, method="arellano"))
| Variable dépendante : dvexptot | ||
| Modèle EF | Modèle EA | |
| age | -6.064 | 4.104*** |
| (7.195) | (0.284) | |
| sex2 | -39.470** | 62.651*** |
| (15.847) | (7.476) | |
| educyr | 26.313*** | |
| (1.262) | ||
| inctot | 0.001*** | 0.002*** |
| (0.0003) | (0.0001) | |
| hinotcov2 | -37.540 | -97.005*** |
| (27.801) | (7.553) | |
| workev2 | -19.919 | 6.564 |
| (33.522) | (12.008) | |
| himcare2 | 242.199*** | 74.702*** |
| (80.002) | (14.360) | |
| Constant | -347.061*** | |
| (22.091) | ||
| Note: | p<0.1; p<0.05; p<0.01 | |
Pour toutes les variables fixes dans le temps, nous n’avons donc pas
de coefficient dans le modèle à effets fixes, contrairement à un modèle
à effets aléatoires. De plus, nous constatons que dans le modèle à
effets aléatoires, peu de coefficients sont statistiquement
significatifs, seuls inctot et
himcare sont significatifs au seuil de 1%, alors que dans
notre modèle à effets aléatoires, presque tous les coefficients sont
significatifs. Cependant, prenons la variable himcare, qui
est significative dans les deux modèles, on constate qu’elle a un
écart-type très important dans le modèle à effets fixes (le coefficient
est également bien plus élevé).
| EA | EF | |
| age | 4.104*** | -6.064 |
| (0.284) | (7.195) | |
| sex2 | 62.651*** | -39.470** |
| (7.476) | (15.847) | |
| educyr | 26.313*** | |
| (1.262) | ||
| inctot | 0.002*** | 0.001*** |
| (0.0001) | (0.0003) | |
| hinotcov2 | -97.005*** | -37.540 |
| (7.553) | (27.801) | |
| workev2 | 6.564 | -19.919 |
| (12.008) | (33.522) | |
| himcare2 | 74.702*** | 242.199*** |
| (14.360) | (80.002) | |
| Constant | -347.061*** | |
| (22.091) | ||
| Note: | p<0.1; p<0.05; p<0.01 | |
Il n’y a donc apriori aucun doute sur le modèle à choisir, c’est le modèle à effets aléatoires car les coefficients sont significatifs, les écarts-types sont plus faibles et, en plus, nous avons des coefficients pour les variables fixes dans le temps. Cependant il ne faut pas oublier que pour un modèle à effets aléatoires, nous devons respecter cette hypothése \(Cov(c_i,x_it) = 0\) si nous voulons un modèle convergent. L’avantage d’un modèle à effets fixes c’est qu’il est convergent si \(Cov(c_i,x_it) ≠ 0\), et ici on peut penser que cette corrélation n’est pas nulle (corrélaltion entre l’éducation et le revenu notamment). Afin d’être sûr, nous pouvons réaliser un test d’Hausman, cela nous permettra de nous aider à faire notre choix entre le modèle à effets fixes ou à effets aléatoires.
\[\begin{cases} H_{0}: Cov(c_i,x_{it}) = 0\\ H_{1}: Cov(c_i,x_{it}) ≠ 0 \end{cases}\]
library(lmtest)
test<-phtest(pool_ea,fixed)
tabtest <- rbind(c("Statistique","$p-value$","$ddl$","Alternative"),
c(round(test$statistic[[1]],2), round(test$p.value,5),
test$parameter[[1]], test$alternative))
tabtest %>%
kbl(caption = "Test d'Hausman", booktabs = T, linesep = "", ) %>%
kable_minimal(full_width = T)
| Statistique | \(p-value\) | \(ddl\) | Alternative |
| 28.77 | 7e-05 | 6 | one model is inconsistent |
On rejette donc \(H_O\) car la p-value est inférieur à 5%. Donc selon le test d’Hausman \(Cov(c_i,x_it) ≠ 0\). On va donc préferer un modèle à effets fixes afin d’avoir un modèle convergent.
On va alors préférer choisir le modèle à effets fixes. En effet, selon le test réaliser précedemment il y aurait une corrélation entre les effets individuels et nos variables observés.
| Dependent variable: | |||
| poolé | EA | EF | |
| age | 4.159*** | 4.104*** | -6.064 |
| (0.283) | (0.284) | (7.195) | |
| sex2 | 63.037*** | 62.651*** | -39.470** |
| (7.492) | (7.476) | (15.847) | |
| educyr | 26.216*** | 26.313*** | |
| (1.264) | (1.262) | ||
| inctot | 0.002*** | 0.002*** | 0.001*** |
| (0.0001) | (0.0001) | (0.0003) | |
| hinotcov2 | -97.688*** | -97.005*** | -37.540 |
| (7.426) | (7.553) | (27.801) | |
| workev2 | 6.474 | 6.564 | -19.919 |
| (11.990) | (12.008) | (33.522) | |
| himcare2 | 72.639*** | 74.702*** | 242.199*** |
| (14.296) | (14.360) | (80.002) | |
| Constant | -349.614*** | -347.061*** | |
| (22.070) | (22.091) | ||
| Note: | p<0.1; p<0.05; p<0.01 | ||
On remarque que dans le modèle à effets fixes, il y a encore la
variable sex, cela signifie alors que certains individus
changent de sexe au cours du temps. Cela semble venir
d’une erreur de remplissage de la base de données, nous pouvons retirer
ces individus et observer s’il y a un changement sur le modèle.
progressite <- data %>%
group_by(id) %>%
summarize(n_unique_sex = n_distinct(sex)) %>% filter(n_unique_sex > 1)
data$SexC <- ifelse(data$id == 21609 | data$id == 45492 | data$id == 48018, "Change","NoChange")
data <- data[data$SexC == "NoChange",]
data <- data[,-11]
fixed2 <- plm(dvexptot ~ age + sex + educyr + inctot + hinotcov + workev + himcare, data = data,
model="within", index = c("id", "year"))
fixed_cor2 <- coeftest(fixed2, vcov=vcovHC(fixed, method="arellano"))
| Variable dépendante : dvexptot | ||
| Modèle EF | Modèle EF corrigé | |
| age | -6.063 | -6.063 |
| (7.167) | (7.195) | |
| inctot | 0.001*** | 0.001*** |
| (0.0002) | (0.0003) | |
| hinotcov2 | -37.533 | -37.533 |
| (33.388) | (27.801) | |
| workev2 | -19.920 | -19.920 |
| (71.695) | (33.522) | |
| himcare2 | 242.199*** | 242.199*** |
| (62.543) | (80.002) | |
| Observations | 103,114 | |
| R2 | 0.001 | |
| R2 ajusté | -1.369 | |
| F Statistic | 6.435*** (df = 5; 43502) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
Hormis le fait que le coefficient associé à la variable
sex n’est plus présent, nous n’observons pas de changement
particulier que ce soit sur les signes des coefficients, leur valeur ou
bien même les écart-types. Désormais les variables qui expliquent les
dépenses dentaires sont inctot &
himcare.
\[Y_t=0.5Y_{t-1}+0.6\epsilon_{t-1}+\epsilon_t \text{ avec } T = 200\]
set.seed(1)
AR1 <-arima.sim(list(ar= 0.5, ma=0.6), n=200)
ggplotly(ggplot(data.frame(AR1), aes(x = 1:200, y = AR1)) +
geom_line(color = "#577559") +
ggtitle("Modèle ARIMA(1,0,1)") +
xlab("Temps") +
ylab("Valeur") +
theme_bw())
acf2y(AR1, 10) %>%
kbl( booktabs = T, linesep = "") %>%
kable_minimal(full_width = T)
| LAG | ACF1 | PACF |
|---|---|---|
| 1 | 0.7155812 | 0.7155812 |
| 2 | 0.3465992 | -0.3390911 |
| 3 | 0.1800755 | 0.2098042 |
| 4 | 0.1199695 | -0.0685944 |
| 5 | 0.1038531 | 0.0838758 |
| 6 | 0.0813634 | -0.0506896 |
| 7 | 0.0536202 | 0.0294945 |
| 8 | 0.0319742 | -0.0203796 |
| 9 | 0.0482315 | 0.0833410 |
| 10 | -0.0022429 | -0.1973629 |
Pour l’ACF, on peut observer un pic à 0.7 environ puis une décroissance, cela veut dire qu’on a une corrélation significative entre les observations, donc le temps présent est corrélé à la dernière période passée ce qui peut indiquer la présence d’un modèle AR dans notre série chronologique.
Pour l’autocorrélation partielle, on a un pic positif à 0.7 ce qui peut indiquer qu’il y a une forte corrélation entre le temps \(t\) et \(t-1\), et donc un modèle autorégressif AR serait adapté. On observe aussi un pic négatif significatif à 0.3 environ, on peut donc penser qu’il y a une relation négative entre la période \(t\) et la période \(t-2\), ce qui indiquerait la présence d’un modèle à moyenne mobile (MA).
On peut donc conclure qu’il y a probablement une combinaison de modèles AR et MA dans notre série temporelle. Les pics positifs significatifs dans l’ACF indiquent une forte corrélation entre les observations passées et les observations actuelles, ce qui suggère un modèle AR. De même, le pic négatif significatif dans la PACF suggère une relation négative entre les observations actuelles et les observations passées, ce qui mène à un modèle MA. Par conséquent, le modèle ARMA pourrait être un bon choix pour modéliser cette série temporelle.
\[Y_t=0.18 +0.47Y_{t-1}+0.61\epsilon_{t-1}+\epsilon_t\]
ARMA <- arima(AR1,order = c(1,0,1))
| ARMA | |
| ar1 | 0.469*** |
| (0.080) | |
| ma1 | 0.608*** |
| (0.083) | |
| intercept | 0.178 |
| (0.203) | |
| Observations | 200 |
| Log Likelihood | -275.272 |
| sigma2 | 0.913 |
| Akaike Inf. Crit. | 558.544 |
| Note: | p<0.1; p<0.05; p<0.01 |
\[Y_t=0.20 + 0.72Y_{t-1}+\epsilon_t\]
AR <- arima(AR1,order = c(1,0,0))
| AR | |
| ar1 | 0.719*** |
| (0.049) | |
| intercept | 0.197 |
| (0.261) | |
| Observations | 200 |
| Log Likelihood | -293.698 |
| sigma2 | 1.100 |
| Akaike Inf. Crit. | 593.396 |
| Note: | p<0.1; p<0.05; p<0.01 |
\[Y_t=0.17 + 0.88\epsilon_{t-1}+\epsilon_t\]
MA <- arima(AR1,order = c(0,0,1))
| MA | |
| ma1 | 0.882*** |
| (0.039) | |
| intercept | 0.166 |
| (0.135) | |
| Observations | 200 |
| Log Likelihood | -287.852 |
| sigma2 | 1.034 |
| Akaike Inf. Crit. | 581.704 |
| Note: | p<0.1; p<0.05; p<0.01 |
On va d’abord effectuer un test de Box-Pierce pour vérifier si les résidus sont des bruits blancs donc qu’ils soient indépendants et pas corrélés.
\[\begin{cases} H_{0}: \epsilon_t \text{ sont indépendants (pas d'autocorrélation dans les résidus)}\\ H_{1}: \epsilon_t \text{ ne sont pas indépendants (autocorrélation dans les résidus)} \end{cases}\]
Si \(H_0\) est rejetée (p-value < \(\alpha\)), les résidus ne sont pas considérés comme des bruits blancs. Si on accepte \(H_0\) (p-value \(\geq \alpha\)), alors les résidus sont des bruits blancs.
Res <- ARMA$residuals
arm<-Box.test(Res)
ResAR <- AR$residuals
a<-Box.test(ResAR)
ResMA <- MA$residuals
m<-Box.test(ResMA)
tab<-rbind(c("ARMA", round(arm$p.value, 4)),
c("AR", round(a$p.value, 4)),
c("MA", round(m$p.value,4)))
colnames(tab) <- c("Modèle", "p-value")
tab %>%
kbl(caption = "Comparaison", booktabs = T, linesep = "") %>%
kable_minimal(full_width = T)
| Modèle | p-value |
|---|---|
| ARMA | 0.8616 |
| AR | 5e-04 |
| MA | 8e-04 |
Les résidus du modèle ARMA(1,1) ont des p-value \(> 0.05\) donc c’est un bruit blanc. Etant donné que les p-value des modèles AR et MA sont inférieurs à \(0.05\), on en conclut que leurs résidus ne sont pas des bruits blancs.
On vérifie que les résidus proviennent d’une distribution normale avec un test de normalité (ici le Shapiro test). Si les résidus suivent une distribution normale alors les erreurs dans les modèles sont réparties de manière aléatoire et ne sont influencées par aucune tendance. Cela assure que les conclusions basées sur le modèle sont fiables et robustes.
\[\begin{cases} H_{0}: \text{les résidus suivent une distribution normale}\\ H_{1}: \text{les résidus ne suivent pas une distribution normale} \end{cases}\]
Si la p-value \(\geq 0.05\), alors on accepte \(H_0\) et les résidus suivent une distribution normale.
arm_shapiro <- shapiro.test(Res)
ar_shapiro <- shapiro.test(ResAR)
ma_shapiro <- shapiro.test(ResMA)
tab2 <- rbind(c("ARMA", round(arm_shapiro$p.value, 4)),
c("AR", round(ar_shapiro$p.value, 4)),
c("MA", round(ma_shapiro$p.value, 4)))
colnames(tab2) <- c("Modèle", "p-value")
tab2 %>%
kbl(caption = "Comparaison", booktabs = T, linesep = "") %>%
kable_minimal(full_width = T)
| Modèle | p-value |
|---|---|
| ARMA | 0.4267 |
| AR | 0.8647 |
| MA | 0.1195 |
Ici, tous les modèles ont des résidus provenant d’une distribution normale.
L’autre critère qu’on pourrait utiliser est le critère d’informations d’Akaike (AIC).
Ce critère mesure la qualité de prédiction d’un modèle en comparant son erreur de prédiction aux informations apportées par son nombre de paramètres.
\[AIC = ln(\frac{SCR_{\epsilon}}{T}) + \frac{2(p+q)}{T}\]
tab3 <- rbind(c("ARMA", round(AIC(ARMA),2)),
c("AR", round(AIC(AR),2)),
c("MA", round(AIC(MA),2)))
colnames(tab3) <- c("Modèle", "AIC")
tab3 %>%
kbl(caption = "Comparaison", booktabs = T, linesep = "") %>%
kable_minimal(full_width = T)
| Modèle | AIC |
|---|---|
| ARMA | 558.54 |
| AR | 593.4 |
| MA | 581.7 |
Plus l’AIC est faible et plus le modèle est de bonne qualité.
Le modèle qui a l’AIC le plus faible est considéré comme étant plus proche de la réalité, ici c’est le modèle ARMA.
Le modèle choisit est cohérent car à la question 1 on utilise aussi un modèle ARMA qui est un modèle à la fois autorégressif et à moyenne mobile. Le fait que l’AIC montre que le meilleur modèle soit l’ARMA(1,1) vient renforcer la pertinence de ce modèle, car il prend aussi en compte les relations dans les données passées et les erreurs.
Le modèle ARMA est le plus adapté pour représenter notre simulation réalisée à la question 1, ce qui est logique étant donné que la simulation a été effectuée à partir d’un modèle ARMA.
Comment synthétiseriez vous vos résultats concernant les dépenses en soin dentaire ?
Avec la décision prise à la question précédemment on peut dire que les dépenses en soin dentaire s’explique par trois variables, le
sex, la possesion ou non de l’aide Medicaire et du revenu des individus. Plus un individu est riche et plus il pourra allouer d’argent à ses dépenses dentaires. L’aide medicare permet également d’augmenter les dépenses en soin dentaire.Cependant on peut-être interloqué de la présence du sexe comme variable explicative alors que nous avons décider de prendre un modèle à effets fixes.